/*!
 **************************************************************************
 * \file convdiffLES.cpp
 * This code solves turbulent natural convection in 3D.
 ***************************************************************************
 * \author    Luis M. de la Cruz [ Sat Aug 23 12:06:36 CDT 2009 ]
 ***************************************************************************
 * The equations to be solved are of the form:
 * \f[
 * \frac{\partial \phi}{\partial t} + \frac{\partial}{\partial x_j} 
 * \Big(u_j\phi\Big) =
 * \frac{\partial}{\partial x_j}\left(\Gamma\frac{\partial\phi}{\partial x_j}\right) + S,
 * \,\,\,\,\,\,\,\,\textrm{for  } j = 1,2,3,
 * \f]
 *
 * where \f$ \phi \f$ is a scalar variable \f$T, \rho, u_1, u_2\f$ or \f$u_3 \f$.
 * \f$ \Gamma \f$ is a diffusion coefficient and \f$ S \f$ is the source term.
 *
 * In this case we have an energy equation in terms of temperature \f$T\f$, 
 * coupled with the Navier-Stokes equations. 
 * The coupled equations are solved with the SIMPLEC strategy.
 * Turbulence is solved using Large-Eddy Simulation (LES) and a selective structuture
 * function to account for turbulent viscosity, as described in Lesieur and Metais 1996.
 * 
 * These equations are solved in the prism defined by
 * \f$ x, y \in [0,1] \times [0,1] \times [0,0.5]\f$. 
 * Two opposite walls have constant temperatures but different. 
 * The other walls are adiabatic. The no-slip condition is applied on all walls.
 * Next figure shows the domain and the oscillating boundary conditions imposed.
 *
 * \image html geomturbo.png "Domain and boundary conditions" width=7cm
 * \image latex geomturbo.eps "Domain and boundary conditions" width=7cm
 *
 *  A resume of the results are shown in the next figure:
 * \image html convturbu.png "Turbulent viscosity and temperature. A comparison with experimental results." width=7cm
 * \image latex convturbu.eps "Turbulent viscosity and temperature. A comparison with experimental results." width=7cm
 *
 *  To compile and run this example type the next commands:
 * \par
   \verbatim
   % make
   % ./convdiffLES \endverbatim
 * 
 ***************************************************************************
 */

#include <iostream>
#include <cmath>


#include "Meshes/Uniform.hpp"
#include "Equations/ScalarEquation.hpp"
#include "Equations/Momentum_XCoDi.hpp"
#include "Equations/Momentum_YCoDi.hpp"
#include "Equations/Momentum_ZCoDi.hpp"
#include "Equations/PressureCorrection.hpp"
#include "Schemes/Simplec.hpp"

// LES Schemes
#include "LES/Schemes/CDS_LES.hpp"
#include "LES/Schemes/CDS_XLES.hpp"
#include "LES/Schemes/CDS_YLES.hpp"
#include "LES/Schemes/CDS_ZLES.hpp"
#include "LES/les_ssf.hpp"

#include "Solvers/TDMA.hpp"
#include "Utils/inout.hpp" 
#include "Utils/num_utils.hpp"
#include "Utils/diff_operators.hpp"
#include "Utils/GNUplot.hpp"

using namespace Tuna;

// Scalar field in one dimension
typedef TunaArray<double,3>::huge ScalarField3D;

template<class T_energy, class T_xmom, class T_ymom, class T_zmom, 
	 class T_press, class T_field>
double SIMPLEC(T_energy &, T_xmom &, T_ymom &, T_zmom &, 
	       T_press &, int, double,
	       T_field &, T_field &, T_field &, T_field &,
	       double, double, double );

double errorT, errorX, errorY, errorZ;

int main() 
{
    double left_wall, right_wall;
    double longitud_x, longitud_y, longitud_z;
    int num_nodos_x, num_nodos_y, num_nodos_z;
    double tolerancia, dt;    
    int max_iter,iteraciones_max;
    double Prandtl, Rayleigh;    
    int frecuencia, init;

// Reading data from the input file
    std::ifstream input_cin ("input");
    input_cin >> longitud_x
	      >> longitud_y
	      >> longitud_z
	      >> num_nodos_x
	      >> num_nodos_y
	      >> num_nodos_z
	      >> left_wall
	      >> right_wall
	      >> dt
	      >> iteraciones_max
	      >> tolerancia
	      >> max_iter
	      >> Prandtl
	      >> Rayleigh
	      >> frecuencia;
//
//  Linear System 
//
    SparseMatrix< Diagonal<double, 3> > A(num_nodos_x, num_nodos_y, num_nodos_z); 
    ScalarField3D             b(num_nodos_x, num_nodos_y, num_nodos_z); 
//
//  Mesh and Scalar field definitions
//
    StructuredMesh<Uniform<double, 3> > malla(longitud_x, num_nodos_x,
					      longitud_y, num_nodos_y,
					      longitud_z, num_nodos_z);

    double dx = malla.getDelta(X);
    double dy = malla.getDelta(Y);
    double dz = malla.getDelta(Z);
    int num_vols_x = malla.getExtentVolumes(X);
    int num_vols_y = malla.getExtentVolumes(Y);
    int num_vols_z = malla.getExtentVolumes(Z);

    malla.print();

    ScalarField3D T(malla.getExtentVolumes());
    ScalarField3D p(malla.getExtentVolumes());
    ScalarField3D pp(malla.getExtentVolumes());
    ScalarField3D u(malla.getExtentVolumes());      // u-velocity
    ScalarField3D v(malla.getExtentVolumes());      // v-velocity
    ScalarField3D w(malla.getExtentVolumes());      // w-velocity
    ScalarField3D us(num_nodos_x, num_vols_y , num_vols_z); // u staggered
    ScalarField3D vs(num_vols_x , num_nodos_y, num_vols_z); // v staggered
    ScalarField3D ws(num_vols_x , num_vols_y , num_nodos_z); // w staggered
    ScalarField3D vort1(malla.getExtentVolumes());      // vorticity_x
    ScalarField3D vort2(malla.getExtentVolumes());      // vorticity_y
    ScalarField3D vort3(malla.getExtentVolumes());      // vorticity_z
    ScalarField3D nut(malla.getExtentVolumes());        // turbulent viscosity

    ScalarField3D Tn(malla.getExtentNodes());
    ScalarField3D pn(malla.getExtentNodes());
    ScalarField3D un(malla.getExtentNodes());
    ScalarField3D vn(malla.getExtentNodes());
    ScalarField3D wn(malla.getExtentNodes());
    ScalarField3D nutn(malla.getExtentNodes());
//
// Valor de los nodos en las fronteras.
//
    Range all = Range::all();
    T(T.lbound(firstDim), all, all) = left_wall;   // Left
    T(T.ubound(firstDim), all, all) = right_wall;    // Rigth

// Highly conducting walls: top and bottom

    int bi = T.lbound(firstDim) + 1, ei = T.ubound(firstDim) - 1;
    Range K(T.lbound(thirdDim), T.ubound(thirdDim));

    for(int i = bi; i <= ei; ++i) {
      T(i, T.lbound(secondDim),K) = left_wall-(i*dx - dx*0.5);  // bottom
      T(i, T.ubound(secondDim),K) = left_wall-(i*dx - dx*0.5); // top
    }

    /***
    NumUtils::staggerFromNodesX(us, un);
    NumUtils::staggerFromNodesY(vs, vn);
    NumUtils::staggerFromNodesZ(ws, wn);
    ***/

    ScalarEquation<CDS_LES<double,3> > energia(T, A, b, malla.getDeltas());
    energia.setGamma(0.6); // Turbulent Prandtl Number
    energia.setDeltaTime(dt);
    energia.setDirichlet(TOP_WALL);
    energia.setDirichlet(BOTTOM_WALL);
    energia.setDirichlet(LEFT_WALL);
    energia.setDirichlet(RIGHT_WALL);
    energia.setNeumann(FRONT_WALL);
    energia.setNeumann(BACK_WALL);    
    energia.setUvelocity(us);
    energia.setVvelocity(vs);
    energia.setWvelocity(ws);
    energia.print();

    Momentum_XCoDi<CDS_XLES<double, 3> > momento_x(us, A, b, malla.getDeltas());
    momento_x.setGamma(Prandtl);
    momento_x.setDeltaTime(dt);
    momento_x.setDirichlet(LEFT_WALL);
    momento_x.setDirichlet(RIGHT_WALL);
    momento_x.setDirichlet(TOP_WALL);
    momento_x.setDirichlet(BOTTOM_WALL);
    momento_x.setDirichlet(FRONT_WALL);
    momento_x.setDirichlet(BACK_WALL);
    momento_x.setUvelocity(us);
    momento_x.setVvelocity(vs);
    momento_x.setWvelocity(ws);
    momento_x.setPressure(p);
    momento_x.print();


    Momentum_YCoDi<CDS_YLES<double, 3> > momento_y(vs, A, b, malla.getDeltas());
    momento_y.setGamma(Prandtl);
    momento_y.setDeltaTime(dt);
    momento_y.setDirichlet(LEFT_WALL);
    momento_y.setDirichlet(RIGHT_WALL);
    momento_y.setDirichlet(TOP_WALL);
    momento_y.setDirichlet(BOTTOM_WALL);
    momento_y.setDirichlet(FRONT_WALL);
    momento_y.setDirichlet(BACK_WALL);
    momento_y.setUvelocity(us);
    momento_y.setVvelocity(vs);
    momento_y.setWvelocity(ws);
    momento_y.setPressure(p);
    momento_y.setTemperature(T);
    momento_y.setRayleigh(Rayleigh);
    momento_y.print();

    Momentum_ZCoDi<CDS_ZLES<double, 3> > momento_z(ws, A, b, malla.getDeltas());
    momento_z.setGamma(Prandtl);
    momento_z.setDeltaTime(dt);
    momento_z.setDirichlet(LEFT_WALL);
    momento_z.setDirichlet(RIGHT_WALL);
    momento_z.setDirichlet(TOP_WALL);
    momento_z.setDirichlet(BOTTOM_WALL);
    momento_z.setDirichlet(FRONT_WALL);
    momento_z.setDirichlet(BACK_WALL);
    momento_z.setUvelocity(us);
    momento_z.setVvelocity(vs);
    momento_z.setWvelocity(ws);
    momento_z.setPressure(p);
    momento_z.print();


    PressureCorrection<Simplec<double, 3> > presion(pp, A, b, malla.getDeltas());
    presion.applyBounds(1, num_nodos_x-1, 1, num_nodos_y-1, 1, num_nodos_z-1);
    presion.setUvelocity(us);
    presion.setVvelocity(vs);
    presion.setWvelocity(ws);
    presion.set_du(momento_x.get_du());
    presion.set_dv(momento_y.get_dv());
    presion.set_dw(momento_z.get_dw());
    presion.setPressure(p);
    presion.print();    


    NumUtils::interpolateToNodes(Tn, T);
    NumUtils::interpolateToNodes(pn, p);
    NumUtils::interpolateToNodesUs(un, us);
    NumUtils::interpolateToNodesVs(vn, vs);
    NumUtils::interpolateToNodesWs(wn, ws);
    InOut::writeToFile_DX(Tn, 0, "./Data/temp.", dx, dy, dz);
    InOut::writeToFile_DX(pn, 0, "./Data/pres.", dx, dy, dz);
    InOut::writeToFile_DX(un, vn, wn, 0, "./Data/velc.", dx, dy, dz);


    ofstream sorsumfile("./Data/sorsum"), nusseltfile("./Data/nusselt"),
      errorTfile("./Data/errorT"), errorXfile("./Data/errorX"), 
      errorYfile("./Data/errorY"), errorZfile("./Data/errorZ");
    double sorsum, nusselt;
    int iteracion;
    
    cout << "\n ----- INICIAMOS SIMPLEC ----- \n";

    for(iteracion = init; iteracion <= iteraciones_max; ++iteracion) {

      sorsum = SIMPLEC(energia, momento_x, momento_y, momento_z, presion,
		       max_iter, tolerancia, u, v, w, nut, dx, dy, dz);
      
//
// LES : Selective Structure Function
//
      NumUtils::interX(u, us);
      NumUtils::interY(v, vs); 
      NumUtils::interZ(w, ws);
      
      DiffOperators::vorticity(vort1, vort2, vort3, u, v, w, dx, dy, dz);
      LES_SSF::calcTurbulentViscosity(nut, u, v, w, dx, dy, dz);
      LES_SSF::selectTurbulentViscosity(nut, vort1, vort2, vort3, 
					dx, dy, dz);

      if( iteracion % frecuencia == 0 ) {
	NumUtils::interpolateToNodes(Tn, T);
	NumUtils::interpolateToNodes(pn, p);
	NumUtils::interpolateToNodesUs(un, us);
	NumUtils::interpolateToNodesVs(vn, vs);
	NumUtils::interpolateToNodesWs(wn, ws);
	NumUtils::interpolateToNodes(nutn, nut);
	InOut::writeToFile_DX(Tn, iteracion, "./Data/temp.", dx, dy, dz);
	InOut::writeToFile_DX(pn, iteracion, "./Data/pres.", dx, dy, dz);
	InOut::writeToFile_DX(un, vn, wn, iteracion, "./Data/velc.", dx, dy, dz);
	InOut::writeToFile_DX(nutn, iteracion, "./Data/nutn.", dx, dy, dz);
      }

      //	nusselt = NumUtils::calcNusselt(T, dx/2.0, dy, dz, 0);
	
	nusseltfile << iteracion << "\t" << nusselt << endl;
	sorsumfile << iteracion << "\t" << sorsum << endl;
	errorTfile << iteracion << "\t" << errorT << endl;
	errorXfile << iteracion << "\t" << errorX << endl;
	errorYfile << iteracion << "\t" << errorY << endl;
	errorZfile << iteracion << "\t" << errorZ << endl;
	
	cout << "=> Iter = " << iteracion 
	     << " s = " << sorsum
	     << " e_T = " << errorT
	     << " e_X = " << errorX
	     << " e_Y = " << errorY
	     << " e_Z = " << errorZ
	     << " Nu = " << nusselt << endl;
    }

    return 0;
}


template<class T_energy, class T_xmom, class T_ymom, class T_zmom,
	 class T_press, class T_field>
double SIMPLEC(T_energy &energy, T_xmom &xmomentum, 
	       T_ymom &ymomentum, T_zmom &zmomentum,
	       T_press &pressure, int max_iter, double tolerancia,
	       T_field &u,T_field &v, T_field &w, T_field &nut,
	       double dx, double dy, double dz ) 
{
    double sorsum = 10.0, tolerance = 1e-02;
    int counter = 0;

    int pausa;

    while ( (sorsum > tolerance) && (counter < 20) ) {

	energy.calcCoefficients(nut);
	Solver::TDMA3D(energy, tolerancia, max_iter);
	errorT = energy.calcErrorL2();
	energy.update();

	xmomentum.calcCoefficients(nut);
	Solver::TDMA3D(xmomentum, tolerancia, max_iter);
	errorX = xmomentum.calcErrorL2();
	xmomentum.update();

	ymomentum.calcCoefficients(nut);
	Solver::TDMA3D(ymomentum, tolerancia, max_iter);
	errorY = ymomentum.calcErrorL2();
	ymomentum.update();

	zmomentum.calcCoefficients(nut);
	Solver::TDMA3D(zmomentum, tolerancia, max_iter);
	errorZ = zmomentum.calcErrorL2();
	zmomentum.update();

	pressure.calcCoefficients();
	Solver::TDMA3D(pressure, tolerancia, max_iter);
	pressure.correction();
	sorsum = fabs( pressure.calcSorsum() );

	//	cout << " counter = " << counter << endl;

	++counter;
    }

    return sorsum;
    
}




